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Abstract 

Development of design tools to furnish optimal acoustic environments for lightweight 
aircraft demands the ability to simulate the acoustic system on a workstation. In order to 
form an effective mathematical model of the phenomena at hand, we have begun by study- 
ing the propagation of acoustic waves inside closed spherical shells. Using a fully-coupled 
fluid-structure interaction model based upon variational principles, we have written a finite 
element analysis program and are in the process of examining several test cases. Future 
investigations are planned to increase model accuracy by incorporating non-linear and 
viscous effects. 
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Introduction 


An important problem in the design of commercial aircraft is noise reduction in both 
the cabin and the outside environment. The need for a low level of noise both inside and 
outside the airplane must be balanced with the need to maintain a light structure. We 
have undertaken a program of computational modeling of this system in order to develop 
the modeling tools necessary for design optimization of future aircraft. 

A favorite computational scheme for studying problems of this type is the finite element 
method. The finite element method was first used for multi-disciplinary applications in 
the 1960’s, when it was found to be appropriate for issues including heat conduction, 
hydrodynamic lubrication, and eigenvalue solutions of the Helmholtz equation for fluids 
vibrating in closed spaces [1], p. 274. One of the earliest instances of the employment of 
finite elements in structural-acoustic problems was by the automotive industry. In 1982, 
Nefske et al. [2] performed this kind of analysis for passenger car and truck cabins. In this 
work they present representations of the acoustic modes within the compartments as well 
as frequency response computations using modal synthesis. 

Much other research has been done to establish a method of coupling the behavior of 
an enclosed fluid to the dynamics of the surrounding structure. Usually this is simply done 
by requiring that no fluid enters the solid surface. This implies that the surface normal 
component of the relative velocity between the fluid and the structure must equal zero. 
The equations of the solid and the fluid are thus solved seperately, and the matching of the 
normal velocities at the surface forms the coupling. It is done this way because to start 
with the equations of motion for the whole system and simply try to apply a finite element 
method (such as a Galerkin method) results in difficult, asymmetric matrices. However, 
the problem with separate solution of the equations is that certain terms drop out which 
in fact would not drop out if a more complicated model of the interaction between fluid 
and structure was used. These terms are the mass terms, and their inclusion would make 
for more precise mathematical predictions. 

Instead of starting from the equations of motion, we have started from a variational 
formulation of Hamilton’s principle, and included the fluid-structure interaction as a con- 
straint on the system. This formulation is based upon the work of Kock and Olson [3]. 
In our first case we have ignored gravitational and viscous effects. Though this is not an 
improvement upon previous methods, we wanted to use this simplest case as a method of 
comparison. It should produce exactly the same result as solving the equations separately 
and matching velocities. If it does so, the stage is set for the inclusion of additional terms 
into the statement of Hamilton’s principle which will improve the accuracy of the model. 


Variational Formulation 

Hamilton’s principle is a general statement about the conservation of energy of a system 
of particles (see, e.g. Bedford [4]). It states that the value of the integral over time of the 
kinetic minus potential energy of the system takes a minimum value. Using the common 
notation of variational calculus (see e.g. Courant and Hilbert [5]), this can be written as 

fh 

6 (T-U)dt- 0 

Jt o 
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in the case of a conservative system, where T is the kinetic energy and U is the internal 
energy. By using the principle of virtual work, Hamilton’s principle is easily extended to 
apply to nonconservative systems (in which some force adds energy to the system), and 
also to systems where the generalized coordinates are not independent. Physically, this 
corresponds to the case where the quantities in question are related by some equations, 
and those equations are known as constraints. With the addition of constraints C and 
virtual work terms W, Hamilton’s principle becomes 

ft i 

6 / (T-U + W + C)dt = 0 
Jt o 

For the problem at hand, we will use equation 1, with a virtual work term supplied by a 
surface force on the sphere, a constraint that enforces velocity matching on the boundary, 
and constraints that enforce local and global conservation of mass. 

The solid shell contributes kinetic energy, potential energy and a surface force term. 


( 1 ) 


T s = 



dU_ dU 
dt dt 



f i £ C £ dflg 
Jn, 2 


( 2 ) 


Here Q. s is the solid volume, T s , the solid boundary, U, the displacement vector, e, the 
strain vector, C, the material stiffness tensor, and f s is the surface force acting on the 
solid. 

By taking into account the assumption that the shell is thin, we can change the volume 
integrals to surface integrals. This is due to the fact that we assume the displacements are 
linearly distributed throughout the thickness of the shell, so that the normal to the element 
is preserved - it remains normal and does not change in length under the deformation. 
Then the displacement compoments U = (I7i, U 2 , W) can be written as 


U\(a i, 02 ,C) = «i(ai,« 2 ) + C0i(ai,ar 2 ,O) , 
a 2> C) = u 2(<*l7 <* 2 ) + a 2, 0) , 

W(a 1 ,a 2 ,C)-= u>(ai,a 2 ) , 

where ai,02»C are general curvilinear coordinates, u\,u 2 ,w are the components of the 
displacement of a point on the reference surface and ,j3 2 are rotations. The thin shell 
assumptions also allow us to write the rotations in terms of uj, u 2 , and w. We can write the 
strain vector e in terms of the displacements u by using the strain-displacement relations 
in general curvilinear coordinates. These relations will not be reproduced here as they are 
lengthy and can be found in many references, see e.g. Kraus [6], p. 31. We will write them 
as e = Au, where A is a matrix operator. 

We now substitute all of these relations into equations 2 and integrate over the normal 
component £. Then the contributions of the solid become 


T s = 



Us = 



400 



Q is the matrix C after the integration - it is multipied by some combinations of powers 
of h. This h is the thickness of the sphere, and W s , of course, is unchanged. The solid 
functional, then, is 

n « = jC [/ rj < A “ )TQA “ ,n ’ i - / r> IT (if '^) dr ‘ - Jr, “ ' 7 ’ dT ‘ ] *' 

Taking the first variational, and integrating by parts in the T s term gives 


= f (A u 

Jr, 


fQASudTs + 


2 dt 2 


- / fs 

Jr, 


■ fiu dr s = o. 


The fluid contributes the following terms: 




where Sir is the volume of the fluid, and e is the specific internal energy of the fluid. 
The constraints of local and global conservation of mass, as well as the fluid-structure 
interaction , are part of the fluid functional. Hence, this functional is of the form 

Kf = Jt jp/ e “ \pfv-v~ Ai + V - \ 2 Pf } dSlf 

+ Ir A3 ( - 1? + Un ) dTs + X2M f dt 

where A l5 A 2 , A 3 are Lagrange multipliers, u n is the normal displacement of the fluid bound- 
ary, Vfi is the normal velocity of the fluid boundary, and M f is the mass of the fluid. 

Taking the first variation with respect to v and setting each integral to zero gives a 
series of equations. These equations reveal the true meanings of the Lagrange multipliers. 
This has been done in detail by Kock and Olson [3], p. 467. We will reproduce only their 
results. The Lagrange multiplier A! = if, where ip is the scalar potential (that is, v = V<f ). 
In addition, A 3 = pj if, and A 2 is still unknown, so it shall be renamed A. After integration 
by parts, the fluid functional becomes 


= IJiJa, { Pfe ~l Pf {Vif)1 + 


— \pf + pj dSlf +A Mj dt. 

In the linearized case it is assumed that the change in fluid density is only a small 
pertubation from the inital density p 0 - By taking Pf = Pa + P and linearizing about 
the initial domain, that is, Slj =»• fi^ 0 andr^ ^ ^f 0 i the functional simplifies further. 
Following Kock and Olson ([3], p. 485), the expression becomes 

n ' - r u* h £ (* - ©’ + s 

+ Jr Un ^ ’ 
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where cq is the speed of sound. In this case taking the first variation and setting each part 
to zero gives 



Finite Element Analysis 

We implement the finite element approximation using standard isoparametric elements 
(see, e.g., Huebner et al. [1], Ch. 5) to discretize the domain. Since the boundaries of each 
element will, of necessity, be curved, the boundary itself must be approximated. For each 
coordinate we will write an expression of the form 

m 

C><) = • ( 5 ) 

1 = 1 

Here T7,£,C are the natural coordinates of the element. N{ are the interpolation functions, 
in natural coordinates, and m is the number of nodes per element. We approximate 
all of the coordinates in this way. In addition, we use the same interpolation functions 
Ni to approximate the variables in the problem, that is, u, v, u>, ip, and A. Displacement 
coordinates are now (u,v,w) instead of ( u i , u 2 , u> ) to prevent confusion. 

We will use two-dimensional elements to model the sphere, and three-dimensional 
elements for the fluid. The approximations are of the form 

m 

*>.7=1 

m 

1=1 

m 

M<) = !>.«■ (6) 

1=1 

where Njj is the 3 x 3m matrix of interpolation functions 

/JV, 0 0 ... 4 0 0 \ 

N= 0 Ni 0 ... 0 N m 0 

\ 0 0 Ni ... 0 0 N m ) 

and il is the vector 

?7=[u 1 Vi Wi ... Um Vm Wm ] . 
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The choice of Lagrange polynomials as interpolation functions allows the nodal values 
V{,Wi, (pi, A,- to have the physical meaning of being the value of the function at that 
node. In our problem the nodal values axe functions of time. Since we will be forcing 
the sphere sinusoidally, we will use Fourier transforms to transform out time and solve 
the probelm in the frequency domain. This will, of course, have no effect on the spatial 
integrals. 

To substitute the approximations into the variational principles derived, it is necessary 
to address the issue of differing coordinate frames. Ths interpolation functions are in local 
coordinates 77 , £, £ , whereas the variational form is in global coordinates — cartesian x,y,z 
in the case of the fluid, and a generalized curvilinear a!,a 2 for the solid. To compound 
the problem, in our case all of the nodal coordinates and element information is obtained 
from the neutral file output of the commercial finite element code PATRAN. The PATRAN 
output has the coordinates of all nodes listed in a global cartesian system. So it is necessary 
to undertake various coordinate transformations before the integrals in the variational 
expressions can be evaluated. 

The fluid expressions give no difficulty in this area. Equation 5 gives the form of 
the tranformation from cartesian to natural coordinates. It is known from calculus that 
dxdydz = \J\drjd(,dC , where \J\ is the determinant of the Jacobian matrix of partial 
derivatives. Thus the integrals will transform in the following manner: 

[ f(x,y,z)dxdy dz = f f'(r]^,C)\J\drjd(dC , 

Jv JV 


where the primes denotate the trasformed function and volume. Substituting equations 6 
into the fluid equations 4 yields the following expressions. 

T \^£r f ^NiNjlJldrjdidC + Vi f po V Ni VNj \ J\dy d£ d( 

. 4 -^ dt 2 Jv cq j Jv 

i,j=i *• 

[ paNi(t),0Nj(ri,0\A<>nd(}fVj = 0 , 

at Jy c l 0 at Js J 

£ j^oN^OW^di 6u]=0, 


tY-S v -%\^^^Sv% NAAdr]dii< 

';J = 1 

/ PoNi(T],i)\J\dr) d£ 6\j = 0 . 

Js J 
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Note that this integration is over the area or volume of the element in natural coordinates, 
so that the range of each variable is (—1, 1). Following Kock and Olson [3], we will write 
these integrals in shorthand, where subscripts denote vector and matrix components: 


m 

E 

t,j=l 

m 


* ht - ^ cf - cr 


dt 2 y 


y dt 


dt y 


= o, 




tj-i 

m r 

E 

l,J=l L 


= 0 , 


- A i K x + ^ Cf X - u? K Xu 


= 0. 


The solid elements demand Em additional coordinate transformation. This will not 
be written out in detail here. The strain-displacement relations written in generalized 
curvilinear coordinates are first transformed to a local arc length coordinate system. Then 
they are fairly simple to trEmsform into natural coordinates using a Jacobian matrix. The 
final result is that equation 3 written with indicial notation in terms of a finite element 
approximation becomes 


”* T f j 2 r l 

E “i / H ki A pk Q n A, r N rj \J\d v di + / -^N ki N tj \J\d v d( 

tj=] L Js Js - 


- [ fiN tj \J\dvd( 

*/ s 



where A is the matrix of strain-displacement relations in natural coordinates, 
write this as 


m 

E 

»,j=i 


\U 


(Pu- 

av^ _ 

’ y ^ dt 2 y * 


= 0 


We can 


Then, putting Eill of these equations together and taking the Fourier trEmsform of both 
sides, the linear system to solve becomes 



H 

3 

1 

O 

3 


u 


' o c^« T o ' 


r u 

u 2 

o 

9- 

1 

O 



-f* XU) 

c <pu o c ^a t 


V 


l -k Au o -a' a J 


A J 


. 0 C^ A 0 ; 


U J 



o 

o 

3 

s 


' u 


r F ' 

+ 

0 -M^ 0 


V 

= 

0 


> 

o 

o 

o 
✓ 


l A J 


0 , 


where the vector [u p A ] is the vector of nodal values. The next step is to use a complex 
solver suitable to this system of equations to solve for the nodal variables. We will use a 
linear system solver taken from the fortran subroutine library LAPACK. 
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Summary and Conclusions 


We have formuated and begun to implement a finite element analysis program based 
on variational principles for the problem of sound wave propagation in a spherical shell. 
At the time of the writing of this report results were not yet available. The program 
has yet to be completely debugged, and several more weeks of work are necessary before 
accurate, complete solutions can be obtained. This work will continue over the course of 
the next year as we incorporate additional effects into our model in our effort to develop 
an improved computational method for structural acoustics problems. 
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